home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 11 / CU Amiga Magazine's Super CD-ROM 11 (1997)(EMAP Images)(GB)(Track 1 of 3)[!][issue 1997-06].iso / cucd / graphics / mpimage / si / median.c < prev    next >
C/C++ Source or Header  |  1995-06-29  |  9KB  |  281 lines
  1. /*
  2. ** File:        median.c             Copyright (c) Truda Software
  3. ** Author:      Anton Kruger         215 Marengo Rd, #2, 
  4. ** Date:        March 1992           Oxford, IA 52322-9383
  5. ** Revision:    1.0                  March 1992
  6. ** 
  7. ** Description: Contains an implementation of Heckbert's median-
  8. **              cut color quantization algorithm.
  9. **
  10. ** Compilers:  MSC 5.1, 6.0.
  11. **
  12. ** Note:       1) Compile in large memory model.
  13. **             2) Delete the "#define FAST_REMAP" statement below
  14. **                in order to deactivate fast remapping.
  15. */
  16.  
  17. // Various changes MJP
  18.  
  19. #undef FAST_REMAP
  20. #include <stddef.h>              /* for NULL                   */
  21. #include <stdlib.h>              /* for "qsort"                */
  22. #include <math.h>
  23. #ifndef _M68881
  24. #include <mffp.h>
  25. #else
  26. #include <m68881.h>
  27. #endif
  28. #include <float.h>               /* for FLT_MAX, FLT_MIN       */
  29. #define MAXCOLORS 256            /* maximum # of output colors */
  30. #define HSIZE     32768          /* size of image histogram    */
  31. typedef unsigned  char byte;     /* range 0-255                */
  32. typedef unsigned  short word;    /* range 0-65,535             */
  33. typedef unsigned  long dword;    /* range 0-4,294,967,295      */
  34.  
  35. /* Macros for converting between (r,g,b)-colors and 15-bit     */
  36. /* colors follow.                                              */
  37. #define RED(x)     (byte)((((x)&31)<<3)|(((x)&31)>>2))
  38. #define GREEN(x)   (byte)(((((x)>>5)&31)<<3)|((((x)>>5)&31)>>2))
  39. #define BLUE(x)    (byte)(((((x)>>10)&31)<<3)|((((x)>>10)&31)>>2))
  40.  
  41. typedef  struct {       /* structure for a cube in color space */
  42.    word  lower;         /* one corner's index in histogram     */
  43.    word  upper;         /* another corner's index in histogram */
  44.    dword count;         /* cube's histogram count              */
  45.    int   level;         /* cube's level                        */
  46.    byte  rmin,rmax;   
  47.    byte  gmin,gmax;   
  48.    byte  bmin,bmax;   
  49. } cube_t;
  50.  
  51. static cube_t list[MAXCOLORS];   /* list of cubes              */
  52. static int longdim;              /* longest dimension of cube  */
  53.  
  54. static void Shrink(cube_t * Cube, word *HistPtr);
  55. static void InvMap(word * Hist, byte ColMap[][3],word ncubes, word *HistPtr);
  56. static int  compare(const void * a1, const void * a2);
  57.  
  58. extern void SetMax(int max);
  59. extern void SetCur(int cur);
  60.  
  61. word MedianCut(word Hist[],byte ColMap[][3], int maxcubes, word HistPtr[])
  62. {
  63.    /*
  64.    ** Accepts "Hist", a 32,768-element array that contains 15-bit
  65.    ** color counts of the input image. Uses Heckbert's median-cut
  66.    ** algorithm to divide the color space into "maxcubes" cubes,
  67.    ** and returns the centroid (average value) of each cube in
  68.    ** "ColMap". "Hist" is also updated so that it functions as an
  69.    ** inverse color map. MedianCut returns the actual number of
  70.    ** cubes, which may be less that "maxcubes".
  71.    */
  72.  
  73.    byte        lr,lg,lb;
  74.    word        i,median,color;
  75.    dword       count;
  76.    int         k,level,ncubes,splitpos;
  77.    void        *base;
  78.    size_t      num,width;
  79.    cube_t      Cube,CubeA,CubeB;
  80.  
  81.    /*
  82.    ** Create the initial cube, which is the whole RGB-cube.
  83.    */
  84.  
  85.    ncubes = 0;
  86.    Cube.count = 0;
  87.    for (i=0,color=0;i<=HSIZE-1;i++){
  88.       if (Hist[i] != 0){
  89.          HistPtr[color++] = i;
  90.          Cube.count = Cube.count + Hist[i];
  91.       }
  92.    }
  93.    Cube.lower = 0; Cube.upper = color-1;
  94.    Cube.level = 0;
  95.    Shrink(&Cube,HistPtr);
  96.    list[ncubes++] = Cube;
  97.  
  98.    /*
  99.    ** The main loop follows.  Search the list of cubes for the
  100.    ** next cube to split, which is the lowest level cube.  A
  101.    ** special case is when a cube has only one color, so that it
  102.    ** cannot be split.
  103.    */
  104.  
  105.     SetMax(maxcubes);
  106.     SetCur(ncubes);
  107.    while (ncubes < maxcubes){
  108.  
  109.       level = 255; splitpos = -1;
  110.       for (k=0;k<=ncubes-1;k++){
  111.          if (list[k].lower == list[k].upper)  
  112.                   ;                            /* single color */
  113.          else if (list[k].level < level){
  114.             level = list[k].level;
  115.             splitpos = k;
  116.          }
  117.       }
  118.       if (splitpos == -1)            /* no more cubes to split */
  119.          break;
  120.  
  121.       /*
  122.       ** Must split the cube "splitpos" in the list of cubes.
  123.       ** Next find the longest dimension of this cube, and update
  124.       ** the external variable "longdim", which is used by the 
  125.       ** sort routine so that it knows along which axis to sort.
  126.       */
  127.  
  128.       Cube = list[splitpos];
  129.       lr = Cube.rmax - Cube.rmin;
  130.       lg = Cube.gmax - Cube.gmin;
  131.       lb = Cube.bmax - Cube.bmin;
  132.       if (lr >= lg && lr >= lb) longdim = 0;
  133.       if (lg >= lr && lg >= lb) longdim = 1;
  134.       if (lb >= lr && lb >= lg) longdim = 2;
  135.  
  136.       /*
  137.       ** Sort along "longdim". This prepares for the next step,
  138.       ** namely, finding the median. Use standard lib's "qsort".
  139.       */
  140.  
  141.       base = (void *)&HistPtr[Cube.lower];
  142.       num  = (size_t)(Cube.upper - Cube.lower + 1);
  143.       width = (size_t)sizeof(HistPtr[0]);
  144.       qsort(base,num,width,compare);
  145.  
  146.       /*
  147.       ** Find median by scanning through the cube, and computing
  148.       ** a running sum. When the running sum equals half the
  149.       ** total for the cube, the median has been found.
  150.       */
  151.       
  152.       count = 0;
  153.       for (i=Cube.lower;i<=Cube.upper-1;i++){
  154.          if (count >= Cube.count/2) break;
  155.          color = HistPtr[i];
  156.          count = count + Hist[color];
  157.       }
  158.       median = i;
  159.  
  160.       /*
  161.       ** Now split "Cube" at the median. Then add the two new
  162.       ** cubes to the list of cubes.
  163.       */
  164.  
  165.       CubeA = Cube; CubeA.upper = median-1;
  166.       CubeA.count = count;
  167.       CubeA.level = Cube.level + 1;
  168.       Shrink(&CubeA,HistPtr);
  169.       list[splitpos] = CubeA;               /* add in old slot */
  170.  
  171.       CubeB = Cube; CubeB.lower = median; 
  172.       CubeB.count = Cube.count - count;
  173.       CubeB.level = Cube.level + 1;
  174.       Shrink(&CubeB,HistPtr);
  175.       list[ncubes++] = CubeB;               /* add in new slot */
  176.       SetCur(ncubes);
  177.    }
  178.  
  179.    /*
  180.    ** We have enough cubes, or we have split all we can. Now
  181.    ** compute the color map, the inverse color map, and return
  182.    ** the number of colors in the color map.
  183.    */
  184.  
  185.    InvMap(Hist, ColMap,ncubes, HistPtr);
  186.    return((word)ncubes);
  187. }
  188.  
  189.  
  190. void static Shrink(cube_t * Cube, word *HistPtr)
  191. {
  192.    /*
  193.    ** Encloses "Cube" with a tight-fitting cube by updating the
  194.    ** (rmin,gmin,bmin) and (rmax,gmax,bmax) members of "Cube".
  195.    */
  196.  
  197.    byte        r,g,b;
  198.    word        i,color;
  199.  
  200.    Cube->rmin = 255; Cube->rmax = 0;
  201.    Cube->gmin = 255; Cube->gmax = 0;
  202.    Cube->bmin = 255; Cube->bmax = 0;
  203.    for (i=Cube->lower;i<=Cube->upper;i++){
  204.       color = HistPtr[i];
  205.       r = RED(color);
  206.       if (r > Cube->rmax) Cube->rmax = r;
  207.       if (r < Cube->rmin) Cube->rmin = r;
  208.       g = GREEN(color);
  209.       if (g > Cube->gmax) Cube->gmax = g;
  210.       if (g < Cube->gmin) Cube->gmin = g;
  211.       b = BLUE(color);
  212.       if (b > Cube->bmax) Cube->bmax = b;
  213.       if (b < Cube->bmin) Cube->bmin = b;
  214.    }
  215. }
  216.  
  217. void static InvMap(word * Hist, byte ColMap[][3],word ncubes, word *HistPtr)
  218. {
  219.    /*
  220.    ** For each cube in the list of cubes, computes the centroid
  221.    ** (average value) of the colors enclosed by that cube, and
  222.    ** then loads the centroids in the color map. Next loads the
  223.    ** histogram with indices into the color map. A preprocessor
  224.    ** directive "#define FAST_REMAP" controls whether the cube
  225.    ** centroids become the output color for all the colors in a
  226.    ** cube, or whether a "best remap" method is followed.
  227.    */
  228.  
  229.    byte        r,g,b;
  230.    word        i,k,color;
  231.    float       rsum,gsum,bsum;
  232.    cube_t      Cube;
  233.  
  234.    for (k=0;k<=ncubes-1;k++){
  235.       Cube = list[k];
  236.       rsum = gsum = bsum = (float)0.0;
  237.       for (i=Cube.lower;i<=Cube.upper;i++){
  238.          color = HistPtr[i];
  239.          r = RED(color);
  240.          rsum += (float)r*(float)Hist[color];
  241.          g = GREEN(color);
  242.          gsum += (float)g*(float)Hist[color];
  243.          b = BLUE(color);
  244.          bsum += (float)b*(float)Hist[color];
  245.       }
  246.  
  247.       /* Update the color map */
  248.  
  249.       ColMap[k][0] = (byte)(rsum/(float)Cube.count);
  250.       ColMap[k][1] = (byte)(gsum/(float)Cube.count);
  251.       ColMap[k][2] = (byte)(bsum/(float)Cube.count);
  252.    }
  253.    return;
  254. }
  255.  
  256. int static compare(const void * a1, const void * a2)
  257. {
  258.    /*
  259.    ** Called by the sort routine in "MedianCut". Compares two
  260.    ** colors based on the external variable "longdim".
  261.    */
  262.  
  263.    word        color1,color2;
  264.    byte        c1,c2;
  265.  
  266.    color1 = (word)*(word *)a1;
  267.    color2 = (word)*(word *)a2;
  268.    switch (longdim){
  269.       case 0:  
  270.          c1 = RED(color1),  c2 = RED(color2);
  271.          break;
  272.       case 1:
  273.          c1 = GREEN(color1), c2 = GREEN(color2);
  274.          break;
  275.       case 2:
  276.          c1 = BLUE(color2), c2 = BLUE(color2);
  277.          break;
  278.    }
  279.    return ((int)(c1-c2));
  280. }
  281.